First-principles study of the ferroelastic phase transition in CaCl2 
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First-principles density-functional calculations within the local-density approximation and the 
pseudopotential approach are used to study and characterize the ferroelastic phase transition in 
calcium chloride (CaCl2). In accord with experiment, the energy map of CaCl2 has the typical 
features of a pseudoproper ferroelastic with an optical instability as ultimate origin of the phase 
transition. This unstable optic mode is close to a pure rigid unit mode of the framework of chlorine 
atoms and has a negative Griineisen parameter. The ab-initio ground state agrees fairly well with 
the experimental low temperature structure extrapolated at K. The calculated energy map around 
the ground state is interpreted as an extrapolated Landau free-energy and is successfully used to 
explain some of the observed thermal properties. Higher-order anharmonic couplings between the 
strain and the unstable optic mode, proposed in previous literature as important terms to explain the 
soft-phonon temperature behavior, are shown to be irrelevant for this purpose. The LAPW method 
is shown to reproduce the plane-wave results in CaCl2 within the precision of the calculations, 
and is used to analyze the relative stability of different phases in CaCl2 and the chemically similar 
compound SrCl2. 



I. INTRODUCTION 

' In recent years, ab-initio methods based on density- 
functional theory have been shown to have predictive 

\ power in the field of thermal structural phase transitions. 
In particular, first-principles calculations within the 
local-density approximation (LDA) have been able to re- 
produce ferroelectric and other instabiliticsJn jncmvskite 

' oxides such as BaTiOa, PbTiOg or SrTi03.Bll3'M These 
ab-initio calculations correctly predict a low-symmetry 
distorted ground state and the eigenvector of the asso- 

\ elated distorting mode. Moreover, the ab-initio explo- 
ration of the energy around this ground state for the rele- 
vant degrees of freedom allows, in general, a parametriza- 
tion of the obtained energy map as an effective hamilto- 
nian. Thermal effects investigated through Monte Carlo 
or Molecular Dynamics simulations restricted to this ef- 
fective hamiltonian show a fair agreement with experi- 

\ ment. 

However, the field of temperature-driven ferroelastic 
\ phase transitions, where some spontaneous strain compo- 
nent can be identified with the order parameter (proper 
] ferroelastic) or is bilinearly coupled with it (pseudoproper 
ferroelastic), represents a particular challenge for first- 
principles calculations. Density- functional methods are 
known to have accuracy problems when reproducing elas- 
tic properties, while in ferroelastics the elastic energy 
plays a fundamental role in the structural instability. 

We report here the results of an ab-initio study of the 
ferroelastic instability in CaCl2 based on the density- 
functional theory within the LDA. A brief preliminary 
account was presented in Rcf. |^. To our knowledge, this 
is the first study of this type in a ferroelastic temperature- 
driven system, and can be considered a benchmark of the 



capability of these ab-initio methods to reproduce ferroe- 
lastic structural instabilities. The case of CaCl2 is par- 
ticularly simple: having a rutile-type structure, it under- 
goes a tetragonal-orthorhombic transition rWtiich is well 
characterized as pseudoproper ferroelasticQu The order 
parameter is one-dimensional and the eigenvector of the 
unstable optic mode is fully determined by symmetry ar- 
guments. This allows an investigation of the origin of the 
structural instability by considering only a few degrees of 
freedom of the structure. In comparison with the case of 
the ferroelectric perovskites, another simplifying feature 
is the fact that the unstable mode is non-polar, so it is 
not necessary to include long-range dipolar interactions 
in any subsequent modeling of the relevant degrees of 
freedom as a local-mode effective hamiltonian. 

In this work, the existence of an orthorhombic (slightly 
distorted rutile structure) ground-state of the system is 
explored, and the total energy is parametrized as a func- 
tion of the relevant degrees of freedom. Special emphasis 
is put on a comparison of the resulting energy map with 
a temperature dependent Landau potential, and the re- 
sulting experimental behavior. We also investigate the 
importance of other secondary degrees of freedom. 

The paper is organized as follows. In Sec. II we de- 
scribe the computational details of the first-principles 
calculations. In Sec. Ill we present the ground state 
of CaCl2 obtained from our ab-initio calculations and its 
comparison with experiment. Temperature effects are 
discussed in detail in Sec. IV. In Sec. V we demon- 
strate the importance of secondary degrees of freedom 
in the determination of the ground state of the system 
and present the results obtained from a direct ab-initio 
minimization of the orthorhombic structure of CaCl2. 
The effects of external stresses exerted on CaCl2 are 
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also considered. Most of the calculations were performed 
with a pseudopotential method but the results were cross 
checked against the highly accurate all-electron linearized 
augmented plane wave (LAPW) method. The latter 
was mainly used to study the relative energies between 
the experimental CaCl2 structure and the denser-packed 
fluorite-type structure in both CaCl2 and SrCl2. In Sec. 
VI we discuss the competition between the three differ- 
ent structures (fluorite-type, rutile-type and orthorhom- 
bic CaCl2-type). Finally we draw some relevant conclu- 
sions in Sec. VII. 



relation was treated withiiiJhe LDA using the Perdew 
and Wang parametrizationtZl of the Ceperley and Alder 
data. For tetragonal and orthorhombic structures we 
took a 5 X 5 X 8 Monkhorst-Pack shifted k-point sam- 
pling in the BZ to give converged total energies, and 
al2xl2xl2 sampling for cubic structures. Self- 
consistency was reached when the total energy of consec- 
utive iterations changed by less than 0.01 mhartree. The 
forces were computed according to Kohler et alJla and 
relaxation of the free atomic coordinates was completed 
when the atomic forces were less than 0.5 mhartree/bohr. 



II. DETAILS OF THE AB-INITIO 
CALCULATIONS 

Part of the calculations presented in this paper were 
performed using density-functional theory within the 
pseudopotential approach with plane waves (PW). For 
these the exchange-correlation energy was lavaluated 
within the local-density approximation, iLDApIl!] using 
the Perdew apd Zunger parametrizationO of the Ceper- 
ley and Aldeill3 interpolation formula for a homogeneous 
electron gas. First-principles norm-conserving pseudopo- 
tentials for Ca and CI were generaied using the scheme 
proposed by Troullier and MartinsEj The Ca pseudopo- 
tential was calculated with a nonlinear core correction in 
the 3(^0-20^^0.50^^0.15 non-spin-polarized valence configu- 
ration, with cutoff radii (in bohr) Ted = 1-29, Tcs = 2.66, 
and Tcp = 3.09. The CI pseudopotential was calculated 
in the Zs^^'^^■ip^^^^?>dP^^^ non-spin-polarized valence con- 
figuration with cutoff radii res = 1.72, r^p = 1.48, and 
Ted = 2.21. A plane- wave basis set up to a kinetic energy 
cutoff of 25 ryd and a 4 x 4 x 5 Monkhorst-PacktJ shifted 
k-point sampling in the Brillouin zone (BZ) were con- 
sidered. Both the basis set and k-point sampling were 
successfully tested to givCj-eonverged total energies. A 
modified Broyden methodic was used for the simulta- 
neous relaxation of the ions and the unit cell while the 
symmetry was preserved. Relaxation of the free atomic 
coordinates was considered accomplished when atomic 
forces were less than 0.05 mhartrec/bohr, and the unit- 
cell relaxation stopped when the stresses were smaller 
than 1 MPa. 

The full-potential LAPW method as implemented in 
the WIEN97 programlla was used to obtain the results 
presented in Sec. VI. The "mufhn tin" radius was cho- 
sen to be 2.3 bohr for all atoms involved in the present 
calculations. The maximum number of plane waves is 
determined by the parameter i^mtA'niax, which was 8 in 
our calculations (i?int is the radius of the atomic sphere 
and i^max is the largest reciprocal wavevector used in the 
LAPW basis set). This led to about 1000 plane waves 
to describe the valence and semi-core states. Local or- 
bitals were used to treat the Ca 3s, 3p and CI 3s semi- 
core states. The density of the core electrons was com- 
puted in the crystal potential of each iteration in the 
self-consistent cycles (thawed core). Exchange and cor- 



III. THE FERROELASTIC GROUND STATE OF 
CaCla 

CaCl2 undergoes a second order ferroelastic phase 
transition at 491K from a high-temperature tetragonal 
rutile-type structure (space group P4:2/mnm, Z = 2) to 
an orthorhombic one (space group Pnnm,Z = 2). The 
phase transition is therefore equitranslational, with an 
order parameter of Big symmetry (antisymmetric for op- 
erations 4 and rrixy)- This is in fact the symmetry of the 
orthorhombic spontaneous strain e = (ei — e2)/2. The 
tetragonal high-temperature structure has also a single 
optical mode of Big symmetry, which according to spec- 
troscopic observations strongly softens as the temper^, 
ature is lowered and causes the structural instability.Q 
Thus, the case of CaCl2 is a textbook example of a pseu- 
doproper ferroelastic transition, where the actual order 
parameter Q is the amplitude of an optical mode, while 
the spontaneous strain, being of the same symmetry, is 
bilinear ly coupled to it. 

We started the ab-initio analysis by searching for the 
tetragonal rutile-type structure of minimal energy, hence- 
forth called the (ab-initio) parent structure. The atomic 
positions, constrained within the symmetry P42 /mnm 
(a single free atomic parameter), and the ratio c/a of the 
tetragonal unit cell were relaxed for different unit cell 
volumes. For each volume the total energy of the relaxed 
structure was calculated and the corresponding results 
were fitted using the Murnaghan equation of statell^ to 
obtain the volume of the P 4,2 /mnm structure of min- 
imal energy. The calculated parent structure is com- 
pared in Table ^ with the expecimental high-temperature 
structure extrapolated at K.E2I The agreement between 
PW and experiment is remarkable considering that LDA 
is known to underestimate unit cell volumes and ther- 
mal effects influence the experimental high-temperature 
structure. Also shown in the table are the structural 
parameters obtained using the LAPW method. The de- 
viation between PW and LAPW, both using the LDA, is 
most likely caused by the different treatment of the Ca-3p 
semi-core states: indirectly through the use of non-linear 
core corrections in PW, and properly in LAPW. 

A sketch of the structure is shown in Fig. |l|. Each 
calcium cation is surrounded by a distorted octahedron 
of chlorine anions of mmm symmetry. Ideal octahedral 
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TABLE L Computed structural parameters and bulk mod- 
ulus for tetragonal CaCl2 compared with the experimental 
values extrapolated to K (Ref. ^). The x atomic coordi- 
nate of the CI atom at 500 K is given in lattice units. We 
include both PW and LAPW results for comparison. 





PW 


Experiment 


Difr.(%) 


LAPW 


a(A) 


6.302 


6.343 


1 


6.242 


c(A) 


4.117 


4.142 


1 


4.068 


x{C\) 


0.3028 


0.3040 


0.4 


0.3032 


Bo(GPa) 


39 






41 




FIG. 1: Sketch showing the parent structure of CaCl2 and 
the eigenvector of the unstable Big mode viewed down the z 
axis. The mode is close to a rotation around the z axis of the 
CaCle octahedra as rigid units. 

coordination would require c/a = 2 - ^/2 = 0.5858 and 
X = {2 — \/2) /2 — 0.2929 as compared with values of 
0.6533 and 0.3028 for our parent structure. This typical 
distortion of anion octahedra in rutile type structures has 
been explained as a simple mechanism for maximizing 
the anion-anion distance, which minimizes its repulsion 
energy (maximum volume per unit cell), while keeping 
fixed the cation-anion distances.Ell Indeed, the unit cell 
volume of the rutile structure would be maximal under 
these constraints for c/a = 0.6325 and x = 0.3000,tll in 
fair agreement with the values calculated ab-initio here 
and those experimentally observed in CaCl2 and-|0ther 
rutile structures, including rather less ionic ones.E3 

The stability of this ab-initio rutile parent structure 
with respect to zone-center distortions was then investi- 
gated by computing the normal mode frequencies at the 
r point. All r modes have real frequencies except the 
mode of Big symmetry, confirming the latter's intrinsic 
instability (see Table |l|) . The frequencies of the stable 
modes agree fairly well with experimental frequencies; de- 
viations of 10% are to be considered normal if it is taken 
into account that we are comparing K calculations with 
experimental values in a structure stabilized at high tem- 



TABLE II: Calculated frequencies (given in cm~^) of all zone- 
center phonons of tetragonal CaCl2 compared with the avail- 
able (Raman-active) experimental frequencies (taken from 
Ref. |7|). 



irrep Theory Exp.(600K) Diff.(%) 

~Arg 2291 2031) 12 

A2g 92.7 

Big 31.5i 17.1 

B2g 261.2 246.7 6 

Eg 169.8 150.0 13 

A2u 248.7 

Biu 71.0, 257.5 

En 76.6,231.6,271.0 



peratures. The Big mode has rather low frequencies in 
many other rutile-type materials and often exhibits a 
strong temperature variatjcta, but in most systems is sta- 
ble at all temperatures.C^EiJft,same compounds it can 
be destabilized via pressure. Ej'E£IlZI Obviously, in the case 
of CaCl2, the "soft" behavior of this mode is so extreme 
that it is unstable at low temperatures. Our calcula- 
tions, which in principle explore the configuration space 
at K, agree with this experimental fact. The tendency 
of the Big mode to be rather soft can be qualitatively 
explained by its "quasi rigid-unit" nature. Indeed, the 
mode is close to a pure rigid unit rotation of the columns 
(parallel to [001]) of distorted chlorine octahedra form- 
ing the structure. Each such column shares corners with 
four neighboring chains (see Fig. 0)r--^o it can be termed 
a "quasi rigid unit mode" (QRUM) ^ilhis can be easily 
checked using the CRUSH package.E3E2l In fact, the Big 
mode does not correspond exactly to a pure rigid unit 
rotation around the z axis due to the different lengths 
of the rotation radius of the chlorine atoms inJjhe CaCle 
octahedra. Ab-initio studies of Karki et alE^ show an 
equivalent situation for the pressure-induced phase tran- 
sition in stishovite. 

Volume effects on the Big mode frequency were also 
calculated. The mode has a negative Griineisen parame- 
ter of around —15, which is in accord f'^if\ experimental 
values in other rutile-type compounds .HE3o The result- 
ing decrease of the mode frequency with pressure is also 
a typical feature of RUMs and QRUMs, and is the rea- 
son why in other rutile-type compounds like rutile itself 
or stishovite the Big instability can be attained at high 
pressures. 

Fig. |2| shows the total energy as a function of 
the amplitude Q of the mode Big, while maintaining 
the rest of the structural parameters at their parent- 
structure values. The eigenvector of this mode (see 
Fig. |l|) has the form 1/V8(000; 000; 110; 110; 110; 110), 
where the displacement vectors are listed in the 
order (e(Ca^), e(Ca"), 6(01^), e(Cl"), e(Cl"i), 6(01^^)) 
(atomic labels are indicated in Fig. |l|). The resulting 
double-well energy map with a minima of 0.44 mhartree 
can be well fitted using a minimal expansion E = Eq + 
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FIG. 2: Total energy (per unit cell) of CaCb as a function of 
the amplitude Q of the B\g mode. The solid line represents 
the double- well fit indicated in the text. 



l/2K(5^-|-l/4a(5'', with a mode amplitude of 1.15 bohr at 
the well minima. The search for the actual ground state 
requires, however, an exploration of the energy of the sys- 
tem as a function of both Q and the orthorhombic strain, 
with which it can be bilinearly coupled. In general, the 
energy map E{Q, e) around the parent tetragonal config- 
uration can be expanded in the form: 



E — En 



-aQ^ + -Coe2+70e + 7'Q'e. (1) 



Here Cq is the bare elastic constant of the orthorhom- 
bic strain e and corresponds to 2(Cii — C12). We have 
included only harmonic terms for the elastic energy, but 
a higher-order coupling term Q^e is also included as it 
will be shown below to be significant. All coefficients in 
this energy expansion can be determined from ab-initio 
calculations: n and a from the above mentioned fit of 
the double well of Fig. ^ and the bare elastic constant Co 
from the parabolic fit of the total energy calculated as a 
hmction of the strain e at Q = (see Fig. ^). The cou- 
pling coefficients 7 and 7' can be obtained from the fit of 
the calculated e-conjugate stress dE/de (cti — (72) vs. Q 
at e = (see Fig. ||b) to the expected behavior 7Q+7'Q'^. 
In Fig. Hb, deviations from the linear behavior are clearly 
observed, so the higher non-linear coupling Q^e is quite 
significant. Table HI lists the values obtained for these 
various coefficients. The resulting energy map E{Q,e), 
has three stationary points: the saddle point at Qo = 0, 
eo = 0, corresponds to the unstable tetragonal structure, 
and the other two symmetry-equivalent stationary points 
locate the ground state of the structure at the twin re- 
lated absolute minima {Qm,£m) = ±(1.27 bohr, 0.029). 
The depth of the minima with respect to the tetragonal 
parent structure is 0.70 mhartree. When compared with 
the relative minimum along Q, it means that nearly 50% 
of the energy minimization corresponds to strain relax- 
ation. This should be compared with the experimental 
Big distortion extrapolated to K. According to Ref. |2^, 
the strain e at room temperature is 0.014, and the exper- 
imental curve e(T) can be extrapolated to e(0) = 0.022. 
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FIG. 3: Upper panel: Fit of the calculated total energy for 
different strained unit cells to obtain the bare elastic constant 
Cq. Lower panel: Calculated ui —(T2 stress as a function of the 
amplitude Q of the Big mode at e = 0, and the corresponding 
fit to -iQ + iQ'. 



TABLE III: Parameters of the total energy expansions 
S((5,e) given in Eq. (§) and Eq. (§). 







_Bit,-mode 








-1.33 


mh/bohr^ a 


1.01 


mh/bohr'* 






Aig-mode 






s 


70.42 


mh/bohr^ 






7 
a 
a' 
c' 


-16.6 
50 

-1.07 
16.8 


Couplings 
mh/bohr 7' 
mh/bohr 6 
mh/bohr^ h' 
mh/bohr^ 


0.80 
-194 
32.1 


mh/bohr"^ 

mh/bohr 

mh/bohr^ 






Elastic 






Co 
C3 


0.68 
3.15 


h (18GPa) C+ 
h (83GPa) C+3 


6.95 
2.15 


h (184GPa) 
h (57GPa) 



Assuming an analogous behavior for Q(T), as expected in 
pseudoproper ferroelastics, the experimental value of 0.85 
bohr for Q at room temperaturecj can be extrapolated 
to 1.33 bohr at K. This value is in excellent agreement 
with the ab-initio value while, on the other hand, the or- 
thorhombic ab-initio strain is about 30% overestimated. 
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IV. TEMPERATURE EFFECTS 

The calculated energy map E{Q,e) can be used as 
the starting point for the "prediction" of thermal behav- 
ior, including the phase transition. This can be done, 
for instance, through the construction of an effective 
hamiltonian consistent with the energetics described by 
E{Q,e), and its further analysis through Monte Carlo 
simulations.u As a zeroth-order approximation, however, 
we can already perform some sensible comparisons with 
experimental results, if we note that the energy map 
E{Q,e) can be taken as an approximation for the free- 
energy (for given Q and e, in the sense of a Landau po- 
tential) at K, and we make the (rather drastic) as- 
sumption that the temperature effects are restricted to 
the thermal renormalization up to positive values of the 
quadratic k coefficient, so that it follows a Landau type 
law K oc {T — To). In other words, the coefficients in Ta- 
ble m can be considered an estimation of the coefficients 
in the Landau free energy, except for k. For instance, ac- 
cording to Eq. (U) , the actual observable elastic constant 
corresponding to e should soften according to the law: 
C{T) = Co — 7^/K(r), so that the phase transition takes 
place when C(T) becomes zero at k{Tc) = j'^/Cq. Ac- 



cording to the parameter values of Table III , this would 
correspond to a value of k{Tc) — 0.41 mhartree/bohr^ 
which, taking into account the effective mass of the 
mode (toci in our Q units) represents a mode frequency 
LOB-ig — 17 cm^^, in good agreement with the experimen- 
tal valuel] of 14 cm~^. 

An expansion of the free energy including a higher- 
order coupling term Q'^e as that given in Eq. (|l|) for 
E(Q,e), was proposed to study the temperature behav- 
ior of the soft optic ghonon observed by Raman spec- 
troscopy experiments.l3 In particular, this coupling term 
was included to explain the large value of the slope ra- 
tio r of iLJ^{T) below and above the phase transition, as 
compared to the value of —2 expected from the simplest 
Landau expansion. By using this extended free-energy 
expansion, it was demonstrated that the slope ratio of 
a;^(T) is no longer —2 (the mean-field value), but a func- 
tion which depends on the free-energy expansion coeL 
ficients: r — {k — 4Q:)/2(a — fc), where k = i^Y/Coli 
Using this result, the calculated non-linear coupling co- 
efficient 7' predicts a slope ratio of a;^(T) of —1.9, which 
is not only quite far from the observed r of —6.5, but 
it corrects the classical value —2 in the opposite direc- 
tion. Trying to explain such a slope ratio through this 
non-linear term in the free energy would require a co- 
efficient 7' one order of magnitude larger than the ab- 
initio calculated value and of opposite sign, which seems 
difficult to accept even considering any reasonable ther- 
mal renormalization effect. Therefore, our results in- 
dicate that this non-classical slope ratio must have a 
quite different origin. In this respect, it is quite inter- 
esting to compare the values of at K extrapolated 
from the experimental linear behavior of uP'{T), below 
Tc (see Fig. ||), with the theoretical expected value ob- 
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FIG. 4: Sketch of the experimental temperature behavior of 
the square of the frequency of the soft mode in CaCl2 (Ref. |^. 
The asterisk and the black circle at K represent ab-initio 
values of lo'^ for Q=0, e—Q and at the absolute minima of the 
total-energy map, respectively. 



tained from the curvature along Q at the absolute min- 
ima of the ab-initio energy map. The experimental value 
is 3.74 mhartree/bohr^ (2790 cm~^) and the theoreti- 
cal one K -I- SaQ^ -I- 67'emQm = 3.73 mhartree/bohr^ 
(2780 cm~^), a surprising excellent agreement. This sug- 
gests that the strong slope change of the linear tempera- 
ture behavior of the mode frequency has a rather simple 
origin: the Q-stiffncss coefficient approaches linearly the 
ground state value once the system is below Tc. This im- 
plies, in general, a correction to the usual Landau descrip- 
tion through a change of slope of the linear decrease of the 
quadratic order parameter coefficient around Tc. This 
kind of behavior has been actually observed in Monte 
Carlo simulations of the model wh|eH|-dfiSjCribed by 
means of an effective Landau potential.Eal2a'EZl The cur- 
vature along Q at the tetragonal maximum is —990 cm^^, 
indicated in Fig. ^ with an asterisk. This is four times 
lower than the extrapolated value from the experimen- 
tal linear behavior of w^(r) above Tc. The deviation of 
this extrapolated value from the actual K value is also 
normal in simplified models like the and is related to 
the displacive degree of the transition. Only in the- 
placive limit would one expect both values to agree. 
A deviation factor of four would correspond to a quite 
normal and realistic intermediate regime where couplings 
and ground state energy scales are of the same order of 
magnitude. 
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V. SECONDARY DEGREES OF FREEDOM 

Until now, in our analysis of the energy we have only 
considered the primary degrees of freedom which deter- 
mine the symmetry of the crystal below the transition, 
i.e., we have expanded the energy only in terms of the am- 
plitude Q of the optic mode Big and the orthorhombic 
strain e. However, we must also consider other secondary 
degrees of freedom which do not determine the change of 
symmetry but are compatible with it and may be impor- 
tant in the study of the ground state of CaCl2. In our par- 
ticular case, these are easily identified as the amplitude 
R of the optic mode Aig and the tetragonal strain com- 
ponents e+ = (ei -I- e2)/2 and 63. The eigenvector of the 
mode Aig has the form 1/^8(000; 000; 110; 110; 110; 110) 
(in the same basis used above for the Big mode). Using 
these new variables, we can write down a new expansion 
for E: 

E = Eo + i^Q' + ^aQ^ + icoe' + jQe + j'Qh + 

+ ISR' + \c+e\ + iCse^ + C+ae+eg + 
+ai?e+ + bRe3 + +a'RQ'^ + b'e+Q^ + c'e^Q^ . (2) 

Our calculations show that allowed fourth order terms 
such as i?^, e*, e^^, e|, Qe"^ are not significant and thus 
have not been included in the previous expansion. We 
have also checked that the contribution of other higher- 
order coupling terms is several orders of magnitude 
smaller than those in Eq. (||). C+, C3 and C+3 are bare 
elastic constants and correspond to 2(Cii -I-C12), C33 and 
2Ci3, respectively. All the coefficients in the expansion 
given by Eq. have been determined from ab-initio cal- 
culations, d comes directly from the phonon frequency 
given in Table ||. The bare elastic constants C+ and C3 
were obtained from the parabolic fit of the total energy 
for Q = versus e+ and £3 at Q = respectively. The 
coefficients a and b result from the calculation of the R- 
conjugate force Fn as a function of e+ and £3 respectively, 
and the elastic constant C+3 from the calculation of ct+ 
(cti +(T2) versus £3. The coupling coefficient a' is obtained 
from the parabolic fit of the force as a function of the 
amplitude Q of the mode Big . b' and c' are obtained by 
fitting Cti and versus the amplitude Q to a parabola. 
Table n Ifl gives a full list of the values obtained for these 
coefficients. 

The resulting energy map can be analyzed in a way 
similar to that used in the case of Eq. p). This is ev- 
ident if we notice that the energy E in Eq. (|) can be 
rewritten in an equivalent form to that given in Eq. (|l|) 
but with the quartic coefhcient a renormalized to a 
smaller value a' given by 0.656 mhartree/bohr'* when 
one introduces the minimum conditions (dE/dR)o = 0, 
{dE/de+)Q = and {dE/des)^ = into Eq. (|). The 
most important contribution to the renormalization of 
a comes from the £+Q^ coupling term, with almost 



TABLE IV; Structural parameters and bulk modulus of the 
fully relaxed orthorhombic phase of CaCb cornpared with the 
experimental values extrapolated to K (Ref. EOl). The x and 
y atomic coordinates for CI atoms at room temperature are 
given in lattice units. 





Computed 


Experiment 


Difrerence(%) 


a(A) 


6.429 


6.446 


1 


b{k) 


6.054 


6.167 


2 


c{k) 


4.088 


4.137 


1 


x{C\) 


0.347 


0.325 


7 


y(ci) 


0.255 


0.275 


7 


Bo(GPa) 


31 







90% of the total value. The energy map has three sta- 
tionary points: the trivial one Qo = 0, £0 = with 
_B = which corresponds to a maximum, and the 
ground state {Qm,£m) = ±(1.54 bohr, 0.033) with R = 
-0.01 bohr, £+ = -0.010, £3 = -0.004 and the depth of 
the minima AE — —0.98 mhartree. As expected, the 
minima are deeper because of a smaller value of the ef- 
fective quartic coefficient. Furthermore, the ground state 
values of Q and e have increased, worsening the agree- 
ment obtained with experiment when secondary modes 
were neglected. In fact, the error in the unit cell vol- 
ume due to the LDA over-binding might be crucial when 
one considers this issue. The value&,of £+ and £3 can be 
compared to the experimental onesu extrapolated to OK 
(£+ = -0.006 and £3 = -0.001). 

The deviation from the classical value of the ratio of 
the slopes of the soft-mode temperature behavior could 
in principle be related to the secondary strain degrees of 
freedom analyzed here. As shown above in general, these 
additional variables renormalize the quartic coefficient a 
to a lower value a' . When considering secondary spon- 
taneous strains clamped at optical frequencies the slope 
ratio becomes r — ~2a/a'. In the present case, for the 
calculated strain coupling strengths, this could mean a 
correction of r to about —3, still very far from the —6.5 
observed. On the other hand, these secondary modes or 
degrees of freedom also affect the calculation done above 
for the curvature along Q of the ab-initio energy map 
at the absolute minima and the corresponding ab-initio 
ground state soft-mode frequency. This becomes 5.34 
mhartree/bohr^ (3980 cm~^), a huge increase of about 
40% with respect to the value obtained disregarding sec- 
ondary modes (see Fig. ^). This is significantly larger 
than the extrapolated experimental OK value. Therefore 
one can conclude that the comparison with experiment 
worsens systematically when secondary strain deforma- 
tions directly related to changes in the unit cell volume 
are included. 

The results obtained above can be compared with the 
outcome of a direct search for the orthorhombic ground 
state of CaCl2. This comparison also leads to an esti- 
mation of the error in the calculations. We minimized 
the total energy of an orthorhombic unit cell with re- 
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spect to the volume, allowing for the relaxation of the 
ratios h/ a and c/a and the atomic positions, constrained 
within the Pnnm space group (two free atomic parame- 
ters). The calculated structure is compared in table IV 
with the experimental orthorhombic structure extrapo- 
lated to K. The agreement is very good considering the 
well known LDA underestimation (the predicted volume 
of the orthorhombic unit cell is Vort = 1075 bohr'^, to 
be compared with the value extrapolated to OK from ex- 
periment (1111 bohr"^)). Using the calculated structural 
parameters listed in tables | and IV we see that in the 
absolute orthorhombic minima e — 0.030, e+ ~ —0.010 
and £3 = —0.007, which agrees fairly well with the 
calculated strains from Eq. (||). From the calculated 
atomic positions of the chlorine atoms we also obtain 
the amplitude of the Big mode, Q = 1.55 bohr and 
the Aig mode, R — —0.06 bohr. The energy difference 
between the orthorhombic and tetragonal structures is 
—0.92 mhartree, to be compared with the depth of the 
minima, AE — —0.98 mhartree, from Eq. (||). Hence, 
these results are consistent with the ground state that 
we have found considering the energy expansion given by 
Eq. (0) with errors of the order of 10%. 




400 420 440 460 480 500 520 540 560 580 
V/molecule(au) 

FIG. 5: Total-energy curves for CaCl2 as a function of the 
unit cell volume per formula unit (f.u.) for the cubic fluorite- 
type (C), tetragonal rutile-type (T) and orthorhombic CaCb- 
type (O) structures. Asterisks mark the volumes per f.u. at 
the minima. The common tangent of the O and C curves is 
shown to highlight the very small energy difference between 
their minima. 



VI. ENERGY- VOLUME CURVES AND 
RELATIVE STABILITY OF DIFFERENT 
STRUCTURES 

We have found it useful to look at the issue of the fer- 
roelastic instability in CaCl2 from the point of view of 
the analysis of energy-volume curves, traditional in the 
field of first-principles calculations. Apart from the con- 
sideration of the two phases involved in the transition, 
a question which comes to mind is whether there might 
be other structures that could compete for stability with 
the experimentally observed CaCl2 structure. A good 
candidate is the cubic fluorite structure, with the anions 
placed in a tetrahedral environment and an eightfold co- 
ordination of the cations, which is the stable arrangement 
for the chemically similar compound SrCl2. 

The calculations in this section are carried out us- 
ing the the highly accurate all-electron LAPW method, 
which for the structural properties of CaCl2 yields es- 
sentially the same results as the PW code, as has been 
shown in Tab. but is in principle more suitable to an- 
alyze small energy differences between structures. We 
have computed the energy- volume curves for CaCl2 in the 
rutile-type (T), the CaCl2-type orthorhombic (O), and 
the fluorite-type cubic (C) structures, and show them in 
Fig. ||. The asterisks under C, O, and T mark the equilib- 
rium volumes per formula unit (468, 520, and 535 bohr'^, 
respectively). 

To put the case of CaCl2 in perspective, we have per- 
formed a similar analysis for SrCl2. Fig. || shows the 
total-energy curves for the same three structure-types as 
for CaCl2. The minimum of the total-energy curve for 
the fluorite-type cubic structure is located at 531 bohr'^. 




550 

V/molecule(au) 

FIG. 6: Total-energy curves for SrCl2 as a function of the 
unit cell volume per formula unit (f.u.) for the cubic fluorite- 
type (C), tetragonal rutile-type (T) and orthorhombic CaCb- 
type (O) structures. Asterisks mark the volumes per f.u. at 
the minima. 



to be compared with the observed experimental value of 
574 bohr'^ . The difference can be attributed to the use of 
the LDA. 

The CaCl2-type orthorhombic structure is very simi- 
lar to the rutile-type tetragonal structure. So, it was to 
be expected that the energy difference is small both in 
CaCl2 and in SrCl2 (0.48 and 0.54 mhartree/f.u., respec- 
tively). In both cases the energies continuously merge as 
the volume increases, since the Big optic mode becomes 
stable and the distortions disappear. This is the signa- 
ture of a group-subgroup symmetry-breaking transition 
in an energy-volume diagram. SrCl2 in the rutile-type 
structure would have an unstable Big mode similar to 
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CaCl2, including the negative Griineisen parameter. 

In CaCl2, the three structures have very similar en- 
ergies. The fluorite-type cubic structure is more stable 
than the rutile-type tetragonal structure (the energy dif- 
ference is 0.38 mhartree/f.u.). So, the cubic and tetrag- 
onal structures are competitive and only the orthorhom- 
bic distortion reduces the energy by an additional 0.1 
mhartree/f.u. below the cubic one. As shown in Fig. 
a very small pressure would be enough to transform 
CaCl2 into the denser fluorite-type cubic structure. This 
small energy difference has led us to re-do the calcula- 
tion with the generalized gradient approximation (GGA) 
for exchange and correlation.E^I In this case the fluorite- 
type structure is about 8 mhartree/f.u. higher than the 
tetragonal one, significantly increasing the critical pres- 
sure. This- finding is consistent with the work of Zu- 
pan et al.c3 on the stishovite/alpha-quartz phase tran- 
sition in Si02, which showed that the compressed phase 
is too stable within the LDA, whereas the GGA corrects 
for this over-binding of the LDA and thus agrees better 
with experiment. In that paper it was argued that it 
is a general feature of the LDA to favor more compact 
(isotropic) structures, while the GGA enhances charge in- 
homogeneities and thus stabilizes the higher-volume and 
lower-symmetry phase. The present case is another ex- 
ample of this general behavior. 

In SrCl2, the ground state of the system corresponds 
to the fluorite-type cubic structure, in accord with the 
experimental observation. The energy difference with re- 
spect to the orthorhombic structure (within the LDA) 
is 9.88 mhartree/f.u., much larger than in the case of 
CaCl2. This can be understood by considering the non- 
bonding forces between the cations. The corresponding 
potential grows steeply at short distances due to strong 
repulsive forces. For this reason the relative energies are 
sensitive to the first-neighbor Ca-Ca and Sr-Sr distances. 
In SrCl2 the cations in the tetragonal structure are very 
close to each other and thus raise the total energy. In 
CaCl2, however, the cations are not so close and thus the 
corresponding energy difference between the structures is 
relatively small. 



VII. CONCLUSIONS 

Ab-initio calculations within the LDA approximation 
without volume corrections are sufficient to reproduce 
the main features of the ferroelastic instability in CaCl2 . 
Our analysis of the energetics of the crystal confirm a 
pseudoproper ferroelastic mechanism with an unstable 
optic mode as the microscopic origin of the experimen- 
tal phase transition. The elastic instability arises due 
to the coupling between this soft Big optic mode and 
the orthorhombic spontaneous strain e. The existence of 
QRUMs in the rutile-type framework structure has been 
shown to play an essential role as well. The calculated 
orthorhombic ground state agrees fairly well with the ex- 
perimental low temperature structure. We have iden- 
tified and calculated the anharmonic terms which have 
relevance in the determination of the ground state. The 
secondary degrees of freedom have been also considered. 
Our energy-map calculations rule out the relevance of 
higher anharmonic couplings proposed in previous liter- 
ature to explain the ratio of the slopes of the squared 
frequencies. Although the system modifies its ground 
state noticeably when volume/strain corrections are con- 
sidered, the LDA approximation has been sufficient for 
obtaining a realistic picture of the microscopic mecha- 
nism of this temperature driven ferroelastic system. 
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